// ///////////////////////////// MIT License //////////////////////////////////// //
//                                                                                //
// Copyright (c) 2010 David Zsolt Manrique                                        //
//                    david.zsolt.manrique@gmail.com                              //
//                                                                                //
// Permission is hereby granted, free of charge, to any person obtaining a copy   //
// of this software and associated documentation files (the "Software"), to deal  //
// in the Software without restriction, including without limitation the rights   //
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      //
// copies of the Software, and to permit persons to whom the Software is          //
// furnished to do so, subject to the following conditions:                       //
//                                                                                //
// The above copyright notice and this permission notice shall be included in     //
// all copies or substantial portions of the Software.                            //
//                                                                                //
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     //
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       //
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    //
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         //
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  //
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN      //
// THE SOFTWARE.                                                                  //
//                                                                                //
// ////////////////////////////////////////////////////////////////////////////// //

//!
//!  Description: recommended fitting algorithm
//!
//!  Levenberg-Marquardt minimization and curve fitting 
//!  Ref:
//!    Is Levenberg-Marquardt the Most Efficient Optimization Algorithm for
//!    Implementing Bundle Adjustment?
//!    Manolis I.A. Lourakis and Antonis A. Argyros
//!    Institute of Computer Science, Foundation for Research and Technology - Hellas
//!    Vassilika Vouton, P.O. Box 1385, GR 711 10, Heraklion, Crete, GREECE lourakis,argyros@ics.forth.gr


#ifndef INTEGRATE_H
#define INTEGRATE_H

#include "structure_io.h"
#include "inout.h"

struct velocity_verlet_control
{
	real_type DT;
	real_type dt;
	real_type t0;
	real_type t1;
};

bool operator >>= (const std::vector<std::string> & input_lines,velocity_verlet_control & get);

struct velocity_verlet
{
	size_t n;
	real_type u;
	vector_array_type  an,a;
	position_t  rn;
	velocity_t  vn;
	
	const mass_t & mass;
	const potential_list_t & pl;
	const constrain_t & fixed;

	velocity_verlet(const mass_t & mass,const potential_list_t & pl,const constrain_t & fixed) : fixed(fixed),pl(pl),mass(mass),n(mass.size()),an(n),a(n),rn(n),vn(n)
	{
	
	}

	bool operator () (real_type & t,real_type dt,real_type DT,
			position_t & r,
			velocity_t & v,
			vector_array_type & f,
			double & u,
			unitcell_t & uc)
	{
		pl.eval(u,f,r);  //! time dependence
		for(size_t i = 0; i < n; i++) a[i]=f[i]/mass[i];

		double te = t+DT;
		double T = 0.0;
	
		//std::cerr << r;

		while(t < te)
		{
			for(size_t i = 0; i < n; i++)
			if(!fixed[i])
			rn[i] = r[i] + dt*v[i] + 0.5*dt*dt*a[i];
			else
			rn[i] = r[i];
			
			pl.eval(u,f,rn);  // at t+dt
			for(size_t i = 0; i < n; i++)
			an[i]=f[i]/mass[i];
			
			for(size_t i = 0; i < n; i++)
			if(!fixed[i])
			vn[i] = v[i] + 0.5*dt*(a[i]+an[i]);
			else
			vn[i] = v[i];

			t+=dt;
			rn.swap(r);
			an.swap(a);
			vn.swap(v);
		}
		return true;
	}


};
std::ostream & operator << (std::ostream & os,const velocity_verlet & get);


struct steepest_descent_control
{
	double kappa;
	double prec; 
	int kmax;
	int every;
};

struct steepest_descent
{
	double kappa;
	double prec; 
	int kmax;
	int every;
	real_type u;
	steepest_descent(const steepest_descent_control & sdc) : kappa(sdc.kappa),prec(sdc.prec),kmax(sdc.kmax),every(sdc.every)
	{}

	size_t n;
	int k;
	void operator () (std::ostream & os,
			//! TODO unit cell input?
			const vector_array_type & r0,
			const structure_t & st)
	{
		std::vector<Eigen::Vector3d> r = r0;
		k = 0;
		n = r.size();
		std::vector<Eigen::Vector3d>  f(n);
		std::vector<Eigen::Vector3d>  rn = r;
		int kk = 0;
		bool stop = false;
		u = 0.0;
		while (!stop && (k < kmax))
		{
			k++;

			st.potential_list.eval(u,f,r);
			for(int i = 0; i < n; i++)
			if(!st.atom_list.fixed[i])
				rn[i] = r[i] + kappa*f[i];

			double d = 0.0;
			for(int i = 0; i < n; i++)
			{
				double max_diff = (r[i]-rn[i]).cwise().abs().maxCoeff();
				if( max_diff > d ) d = max_diff;
			}
			if(d < prec ) stop = true;
			rn.swap(r);

			
			if(++kk > every)  //!TODO t=t0 print
			{
				kk = 0;
				os << n << ' ';
				os << std::endl;
				os << st.unitcell(0,0) << ' ';
				os << st.unitcell(0,1) << ' ';
				os << st.unitcell(0,2) << ' ';
				os << st.unitcell(1,0) << ' ';
				os << st.unitcell(1,1) << ' ';
				os << st.unitcell(1,2) << ' ';
				os << st.unitcell(2,0) << ' ';
				os << st.unitcell(2,1) << ' ';
				os << st.unitcell(2,2) << ' ';
				os << "potential=" << u << ' ';
				os << std::endl;
				for(std::vector<Eigen::Vector3d>::size_type i = 0; i < n; i++)
				{
					os << st.atom_list.label[i] << "  ";
					os << ios_real << r[i].x();
					os << ios_real << r[i].y();
					os << ios_real << r[i].z();

					os << ios_real << ""; // for testing!!!!!!!!!! ONLY /Todo
					os << ios_real << "";
					os << ios_real << "";

					os << ios_real << f[i].x();
					os << ios_real << f[i].y();
					os << ios_real << f[i].z();
					os << std::endl;
				}
			}

		}
				
	}

};

bool operator >>= (const std::vector<std::string> & input_lines,steepest_descent_control & get);
std::ostream & operator << (std::ostream & os,const steepest_descent & get);



#endif //INTEGRATE_H

